% Value Function Iteration for a menu cost model with:
%
% 1. Constant aggregate demand
% 2. I.i.d. normal growth of the price level (non-zero mean)
% 3. AR(1)-normal idiosyncratic productivity shocks with homoskedastic errors
%
% Inputs:
%    rPi: The profit function
%    ErV: An initial guess for ErV (the expected value function)
%    beta: the discount factor
%    K: The menu cost
%    theta: Elasticity of demand
%    rp: The grid for the real price
%    Proba: The probability transition matrix for a (the idiosyncratic shock)
%    probInfl: The probability transition matrix for the real price
%       conditional on no nominal price change
%    tolerance: The error tolerance level of the iteration loop
%    tloop: If one the program prints the time the loop took
%    yesprint: If yesprint == 1 then the function prints 'progress reports'
%
% Outputs:
%    rV: The real valued value function
%    F: The policy function (log price)
%
% Jon Steinsson and Emi Nakamura, May 2006
%********************************************************

function [rV,F, IF] = VFIteration_hom_infl(rPi,ErV,beta,K,theta,rp,Proba,ProbInfl,...
    tolerance,tloop,yesprint)

rpgridsize = size(ProbInfl,1);
agridsize = size(Proba,1);
irp=[1:1:rpgridsize]';

time1 = cputime;
i = 1;
ErV_new = ErV;
ErV_old = zeros(rpgridsize,agridsize);
Vsubopt = 4*tolerance;
while Vsubopt/2 > tolerance
    ErV_old = ErV_new;
    % Expected value if nominal price is not changed
    ErV_notchange = rPi + beta*ErV_old;
    % Expected value if nominal price is changed
    ErV_change = ones(rpgridsize,1)*max(rPi-(theta-1)*K/theta+beta*ErV_old);
    % Compare the two and choose the higher
    ErV_compare = zeros(2,rpgridsize,agridsize);
    ErV_compare(1,:,:) = ErV_notchange; 
    ErV_compare(2,:,:) = ErV_change;
    ErV_new = squeeze(max(ErV_compare));
    % Take expectations
    ErV_new = ProbInfl*ErV_new*Proba';
    
    Vdifmax = (beta/(1-beta))*max(max(ErV_new - ErV_old));
    Vdifmin = (beta/(1-beta))*min(min(ErV_new - ErV_old));
    Vsubopt = Vdifmax-Vdifmin;
   
    if Vsubopt/2 <= tolerance
        Vadj = (Vdifmax+Vdifmin)/2;
    end
    if yesprint == 1
        if mod(i,5) == 0
            fprintf('.')
        end
        if mod(i,400) == 0
            fprintf('\n')
        end
    end
    i = i + 1;    
end
ErV = ErV_new + Vadj;

if tloop == 1
    fprintf('\n\n')
    fprintf('Time of Loop: %8.3f\n\n',cputime - time1)
end

rV_notchange = rPi + beta*ErV;
[rV_change rp_change] = max(rPi-(theta-1)*K/theta+beta*ErV);
rV_change = ones(rpgridsize,1)*rV_change;
rV_dif = rV_change - rV_notchange;

F_change = (rV_dif>0);
F_notchange = (F_change==0);
F = F_notchange.*(rp*ones(1,agridsize))+F_change.*(ones(rpgridsize,1)*rp(rp_change,1)');
IF = F_notchange.*(irp*ones(1,agridsize)) + F_change.*(ones(rpgridsize,1)*rp_change);
rV = F_notchange.*rV_notchange + F_change.*rV_change;
